!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!!
!! Define a state variable class.
!!
!! \n
!!
!! This module defines a state variable abstract type \c c_stv which
!! can be used as a building block for the implementation of linear
!! and nonlinear solvers, time integrators and so on. The main purpose
!! of this type is decoupling the details of the state variable
!! representation, such as internal fields and arrays with an
!! arbitrary number of dimensions, from the implementation of some
!! general purpose algorithms. Moreover, this allows a single \c c_sty
!! variable to be used in various such algorithms, avoiding copies.
!!
!! The operators defined for \c c_stv essentially are those required
!! by a vector space, optionally including a scalar product (Hilbert
!! space)
!!
!! \note All the subroutines and functions should have no side effects
!! and should be declared \c pure. However, there are many
!! restrictions concerning pure procedures which essentially prohibit
!! using them together with polymorphic objects: see section <em>12.7
!! Pure procedures</em> of the Fortran 2008 standard, and especially
!! C1278a in <em>CORRIGENDUM 1</em>; plus the fact that pure procedure
!! can not make any MPI call. For this reason, none of the procedures
!! defined in this module is pure.
!!
!! For further details, see the comments in \c c_stv.
!<----------------------------------------------------------------------
module mod_state_vars

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_state_vars_constructor, &
   mod_state_vars_destructor,  &
   mod_state_vars_initialized, &
   c_stv

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 !> State variable
 !!
 !! As a general rule, for all the type-bound procedures the passed
 !! argument is the one which is set by the procedure (an exception is
 !! the scalar product, which is implemented as a function).
 !!
 !! We avoid functions in the type fields because typically this would
 !! require frequent allocations/deallocations of the function
 !! results, while the subroutine arguments can be allocated once and
 !! then passed around thanks to the \c inout intent.
 !!
 !! In general, when extending the type \c c_stv, both allocatable and
 !! pointer fields will be added. This creates some ambiguity when
 !! creating and copying \c c_stv objects. We consider the
 !! following interface:
 !! <ul>
 !!  <li> the user is responsible for setting all the fields of any \c
 !!  c_stv object he/she defines; none of the functions introduced
 !!  here is meant for initializing a new object
 !!  <li> general purpose subroutines operating on \c c_stv objects
 !!  which require local variables of the same type should define such
 !!  local variables as allocatable and allocate them with
 !!  \code
 !!   allocate( local_var , source=input_arg )
 !!  \endcode
 !!  so that all the fields are set as required: pointers point to
 !!  the same object and allocatable fields have distinct storage; in
 !!  other words, source allocation is the copy constructor
 !!  <li> two already defined objects can be made equal using the user
 !!  defined type bound procedure <tt>c_stv\%copy</tt>.
 !! </ul>
 !!
 !! Note that for efficiency reasons it might be useful to override
 !! the implementation of the various operators provided here.
 !!
 !! We also need a comment concerning the parallel execution. In
 !! general terms, there are two kinds of distributed data: those for
 !! which shared data are assumed to be copies of the same value and
 !! those for which shared data are supposed to be added together. All
 !! the operations defined in \c c_stv are correct as far as all the
 !! operants are of the same type, which then is also the type of the
 !! result. For this reason, there is nothing in \c c_stv which refers
 !! explicitly to this aspect. However, it's the user's responsibility
 !! to avoid making operations with variables of different kinds, in
 !! which case the result is meaningless.
 type, abstract :: c_stv
 contains
  !> \f$x+=y\f$
  procedure(i_incr), deferred, pass(x) :: incr
  !> \f$x*=r\f$
  procedure(i_tims), deferred, pass(x) :: tims
  !> \f$x+=ry\f$ (requires a temporary)
  procedure,                   pass(x) :: inlt
  !> \f$z=x\f$
  !!
  !! Notice that copy works on two objects which are already defined
  !! (allocated) by copying the fields of \c x into the corresponding
  !! fields of \c z. This subroutine will be (in general) different
  !! from the standard assignment operator, because various fields
  !! which are not supposed to change among objects of the same type
  !! do not need to be copied (typically pointers to grids, bases
  !! etc.). In general, this is the function that one wants to use
  !! when defining algorithms operating over \c c_stv objects. For the
  !! same reason, this function is also not equivalent to source
  !! allocation.
  procedure(i_copy), deferred, pass(z) :: copy
  !> \f$z=x+y\f$
  procedure,                   pass(z) :: add
  !> \f$z=rx\f$
  procedure,                   pass(z) :: mlt
  !> \f$z=x+ry\f$
  procedure,                   pass(z) :: alt
  !> \f$x+=\sum_i r_iy_i\f$
  procedure,                   pass(x) :: inlv
  !> \f$z=\sum_i r_iy_i\f$
  procedure,                   pass(z) :: lcb
  !> scalar product
  !!
  !! There are situations where a scalar product is never required: to
  !! simplify the implementation we provide here a dummy version
  !! of this subroutine.
  !!
  !! \note For parallel execution, the scalar product must be returned
  !! on <em>all</em> the processors.
  procedure,                   pass(x) :: scal
  !> show something, for debug
  procedure,                   pass(x) :: show
  !> Copy constructor (compiler bug)
  !!
  !! The contruct <code>allocate(y,source=x)</code> is still not
  !! completely supported: see DPD200180295 and
  !! http://software.intel.com/en-us/forums/showthread.php?t=103821
  !!
  !! The subroutine source is a temporary workaround where the copy is
  !! constructed by the user.
!  procedure(i_source), deferred, pass(x) :: source
! Workaround of the woraround: see the comments to i_copy for details
  procedure, pass(x) :: source
  procedure(i_source_vect), deferred, pass(x) :: source_vect
 end type c_stv

 abstract interface
  subroutine i_incr(x,y)
   import :: c_stv
   implicit none
   class(c_stv), intent(in)    :: y
   class(c_stv), intent(inout) :: x
  end subroutine i_incr
 end interface

 abstract interface
  subroutine i_tims(x,r)
   import :: wp, c_stv
   implicit none
   real(wp),     intent(in)    :: r
   class(c_stv), intent(inout) :: x
  end subroutine i_tims
 end interface

 abstract interface
  subroutine i_copy(z,x)
   import :: c_stv
   implicit none
   class(c_stv), intent(in)    :: x
   class(c_stv), intent(inout) :: z
  end subroutine i_copy
 end interface

 abstract interface
  subroutine i_source(y,x)
   import :: c_stv
   implicit none
   class(c_stv), intent(in)               :: x
   class(c_stv), allocatable, intent(out) :: y
  end subroutine i_source
 end interface
 abstract interface
  subroutine i_source_vect(y,x,m)
   import :: c_stv
   implicit none
   integer, intent(in) :: m
   class(c_stv), intent(in)  :: x
   class(c_stv), allocatable, intent(out) :: y(:)
  end subroutine i_source_vect
 end interface
 
 ! private members

! Module variables

 ! public members
 logical, protected ::               &
   mod_state_vars_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_state_vars'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_state_vars_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.   ) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_state_vars_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_state_vars_initialized = .true.
 end subroutine mod_state_vars_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_state_vars_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_state_vars_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_state_vars_initialized = .false.
 end subroutine mod_state_vars_destructor

!-----------------------------------------------------------------------

 subroutine inlt(x,r,y)
  real(wp),     intent(in) :: r
  class(c_stv), intent(in) :: y
  class(c_stv), intent(inout) :: x

  class(c_stv), allocatable :: tmp

   ! Compiler bug (ifort Version 12.1.1.256)
   ! The following two allocations should be equivalent, however when
   ! using the first one the input variable y gets modified by the
   ! call to mlt, possibly due to some poiter copy. The bug is
   ! DPD200180295 and is discussed here:
   ! http://software.intel.com/en-us/forums/showthread.php?t=103821
      !allocate(tmp,source=x)
      !allocate(tmp,source=y)
      call y%source(tmp)
   call tmp%mlt( r , y )
   call   x%incr( tmp )
   deallocate(tmp)

 end subroutine inlt
 
!-----------------------------------------------------------------------

 subroutine add(z,x,y)
  class(c_stv), intent(in) :: x, y
  class(c_stv), intent(inout) :: z

   call z%copy( x )
   call z%incr( y )

 end subroutine add
 
!-----------------------------------------------------------------------

 subroutine mlt(z,r,x)
  real(wp), intent(in) :: r
  class(c_stv), intent(in) :: x
  class(c_stv), intent(inout) :: z

   call z%copy( x )
   call z%tims( r )

 end subroutine mlt
 
!-----------------------------------------------------------------------

 subroutine alt(z,x,r,y)
  real(wp), intent(in) :: r
  class(c_stv), intent(in) :: x, y
  class(c_stv), intent(inout) :: z

   call z%mlt( r , y )
   call z%incr( x )

 end subroutine alt
 
!-----------------------------------------------------------------------

 subroutine inlv(x,r,y)
  real(wp), intent(in) :: r(:)
  class(c_stv), intent(in) :: y(:)
  class(c_stv), intent(inout) :: x

  integer :: i

   do i=1,size(r)
     call x%inlt( r(i) , y(i) )
   enddo

 end subroutine inlv

!-----------------------------------------------------------------------

 subroutine lcb(z,r,y)
  real(wp), intent(in) :: r(:)
  class(c_stv), intent(in) :: y(:)
  class(c_stv), intent(inout) :: z

  integer :: i

   call z%mlt( r(1) , y(1) )
   do i=2,size(r)
     call z%inlt( r(i) , y(i) )
   enddo

 end subroutine lcb

!-----------------------------------------------------------------------

 !> Scalar product
 function scal(x,y) result(s)
  class(c_stv), intent(in) :: x, y
  real(wp) :: s
   
  character(len=*), parameter :: &
    this_fun_name = 'scal'
  character(len=*), parameter :: &
    err_msg(4) = (/ &
   "Please, consider that you have to provide a real implementation ",&
   "for this function if you want to define a Hilber space.         ",&
   "This function in mod_state_vars is provided only to simplify the",&
   "implementation for cases where no scalar product is required.   " &
                 /)

   call error(this_fun_name,this_mod_name,err_msg)
 end function scal

!-----------------------------------------------------------------------

 subroutine show(x)
  class(c_stv), intent(in) :: x
   call warning('show',this_mod_name,        &
     'Overload this function if you need it.')
 end subroutine show

!-----------------------------------------------------------------------

! compiler bug: the following subroutine should not be used: use
! source allocation instead (see the comments above in this module).
  subroutine source(y,x)
   class(c_stv), intent(in)               :: x
   class(c_stv), allocatable, intent(out) :: y

   ! This subroutine must be ALWAYS overridden!!
  end subroutine source

end module mod_state_vars

